Combined Sampled Parameters

Let’s combine and organize our sampled parameters in a data frame. We’ll use an equation-specific function to ensure completeness and repeatability for multi-well matching.

# The MHSample function will sample and combine all parameters based on their 
# DefineDistribution classes. Method input will return different output type based on 
# requirement to speed up MCMC calculations in the future.. Ignore for now.?
MHSample<-function(n=10000,qi,b,di,dlim){
  #Sample relevant parameters for MH function: qi, b, di, dlim
  #Inputs of class distribution created prior to function call, containing 
  # distribution information
  # for each parameter
  #Function will output a Monte Carlo matrix with n x m dimensions, where n is sample number
  # and m is the number of parameters needed for MH equation
  
  #Step through each parameter and call sample.any function
  qi.out <- sample.any(n,qi)
  b.out <- sample.any(n,b)
  di.out <- sample.any(n,di)
  dlim.out <- sample.any(n,dlim)
  
  #Combine each returned column into dataframe for output
  sampledParams <- data.frame("qi"=qi.out,"b"=b.out,"di"=di.out,"dlim"=dlim.out)
  
  #Return sampled dataframe
  return(sampledParams)
}

# Pass the previously defined distributions into the MH.sample function
# Note this is repeating our earlier sampling in a cleaner form
MH.sampled.parameters <- MHSample(n = 10000,qi = qi,b = b,di = di,dlim = dlim)

# Return in table, displaying first 10 rows only
kable(MH.sampled.parameters[c(1:10),], align = 'c') %>% 
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive",
                                      full_width = F, position = "center"))

Historical Rate Synthesis

Using the combined parameter dataframe and the modified hyperbolic equation, we can define a function to use each row of the parameter dataframe as a Monte Carlo iteration and substitute it into the arps equation across the time interval for the sample well. The function will create and return a mcIteration object for each iteration. This object is a list with slots for Iteration name, Parameters, and Synth. Some additional slots will be created and populated later as we approach cost and forecasting.

# MH.synthesize, will be defined to carry out rate synthesis across the actual time interval.
MH.synthesize<-function(parameters,q.abd,t){
  #Parse up parameters and pass into MH function call
  fc <-  MH(parameters[["qi"]],parameters[["b"]],parameters[["di"]],parameters[["dlim"]],q.abd,t)
  
  #Simplify parameters into a single string representation
  parametersName<-paste(names(parameters),round(parameters,digits=2),sep = "=",collapse=" ")
  
  #Create mcIteration imitation list to improve performance and ease of combining data
  #Return mcIteration. Remaining list items to be appended by later functions
  return(list("Iteration" = parametersName,"Parameters" = parameters,"Synth"=fc))
}

To properly synthesize production data for the sample well, we need to create a time vector from the peak oil rate to the last producing month.

# Define a time interval of actuals to synthesize across, using 100/05-17-081-16W6/00
uwi<-"100/05-17-081-16W6/00"

# Subset production including only 100/05-17-081-16W6/00
prod <- formated_production$CD[which(formated_production$CD$Unique.Well.ID == 
                                              uwi),]

# Find peak oil month for this uwi
uwi.pm <- matrix[["Oil.Peak.Month"]][which(matrix$eachUWI==uwi)]

# Refine production datatable to only include data after peak oil rate
prod_ap<-prod[-c(1:(uwi.pm-1)),]

# Find months on prod
mop <- sorted_matrix[which(sorted_matrix$eachUWI == uwi),"Months.on.Prod"]

# Create time vector starting at 0 for peak oil month and using cumulative sum for days
synth.time <- cumsum(c(0,head(prod_ap$Days.In.Month,-1)))

# Let's apply the MH.synthesize function in an apply call to execute across each row
# of the parameters dataframe
oil.synth<-apply(MH.sampled.parameters,1,MH.synthesize,q.abd = o.abd, t = synth.time)

# Print first object of oil.synth to show mcIteration object
print(oil.synth[[1]])

## $Iteration
## [1] "qi=637.8 b=1.46 di=0.73 dlim=0.07"
## 
## $Parameters
##          qi           b          di        dlim 
## 637.8041908   1.4575283   0.7250933   0.0700000 
## 
## $Synth
##  [1] 637.8042 492.5824 406.2513 349.5694 310.1796 279.2837 255.6979 235.8914
##  [9] 219.4737 206.8542 194.7880 184.6325 175.4065 167.4857 160.1658 153.5850
## [17] 147.8145 142.3815 137.5680 132.9939 128.7763 125.1160 121.4741 118.1879
## [25] 115.0127

Each mcIteration list, defining a Monte Carlo iteration and historical synthesis, can now be evaluated against the actual data to search for an appropriate model.

Soft L1

Soft L1 loss is a robust, continuous variation of the common Huber Loss that gives less weight to outliers than a L2 loss, making it useful for fitting time series production data. Our Soft L1 definition will return the individual loss of each point and an average cost for all residuals input.

soft_L1<-function(r){
  IndvCost <- 2*((1 + abs(r))^0.5 - 1)
  Cost <- sum(IndvCost,na.rm = TRUE) / length(IndvCost)
  return(list(IndvCost,Cost))
}

To apply the loss equation we’ll define a function to apply the Soft L1 loss to the residuals of each model match. Due to the non-linear nature of production data, the loss function will be applied to the log residuals. The function will populate the IndvCost and Cost slots of each mcIteration object. For this illustration, the function will be defined in the background and we can see the result for the first iteration object.

## $Iteration
## [1] "qi=637.8 b=1.46 di=0.73 dlim=0.07"
## 
## $Parameters
##          qi           b          di        dlim 
## 637.8041908   1.4575283   0.7250933   0.0700000 
## 
## $Synth
##  [1] 637.8042 492.5824 406.2513 349.5694 310.1796 279.2837 255.6979 235.8914
##  [9] 219.4737 206.8542 194.7880 184.6325 175.4065 167.4857 160.1658 153.5850
## [17] 147.8145 142.3815 137.5680 132.9939 128.7763 125.1160 121.4741 118.1879
## [25] 115.0127
## 
## $IndvCost
##  [1] 0.5799130 0.6213472 0.7045805 0.7993099 0.8145477 0.8342588 0.7150210
##  [8] 0.8231342 0.8641477 0.9012344 0.9399487 1.0331270 0.9826236 1.0292325
## [15] 1.0180480 1.0188545 1.0867748 1.0848554 1.0857151 1.0637195 1.1046275
## [22] 1.1677599 1.1636212 1.2077019 1.1906771
## 
## $Cost
## [1] 0.9533912

Using the cost function applied to each iteration, we can minimize the cost to find our bestfit model. Using this bestfit, we can also visualize the cost model results for each production point in the fit.

Using the log residuals and the Soft L1 cost function, let’s find the top 100 models from our initial 10,000 iterations and see how the matches vary from the bestfit. We’ll colour them by average cost to compare the forecasts.

Top 1% of Models Using Soft L1 Loss

We can see that the top 100 forecasts do a good job of covering the range of actual oil rates, with the lowest cost forecasts tightly grouped around the actuals. However, one observation is that the forecasts are not evenly distributed above and below the actuals to appropriately capture uncertainty in the forecasts. Let’s explore another loss function to see if we can improve this behaviour.

PDCA Loss

The Soft L1 loss function does a decent job at finding a best fit model, but is there a better loss function that can be used? Let’s look at a loss function defined in SPE-174784-PA (ADD LINK HERE AS REFERENCE) by Fulford et al. that takes a similar approach, but is designed to capture more uncertainty in the production data. This loss function uses the average and standard deviation of the individual cost values across the fit, comparing them to a baseline fit. This will prove particularily useful when we move into the PDCA methods later in the CAST program.

For now, let’s define this loss function as PDCA Loss and see how it improves the fits from our Monte Carlo.

# As our Soft L1 was defined, pdca_loss will return slots for IndvCost and Cost
pdca_loss<-function(r, Eps_mean_min = 0, Eps_std_min = 0){
  IndvCost <- r^2 # pdca loss simply uses L2 loss for individual cost
  
  Eps_mean <- mean(r,na.rm = TRUE) #Mean of residuals
  Eps_std <- sd(r,na.rm=TRUE) #Std of residuals
  
  Cost <- ((Eps_mean - Eps_mean_min)/.1)^2+((Eps_std - Eps_std_min)/.01)^2
  return(list(IndvCost,Cost))
}

Just as we did with the Soft L1, let’s use the PDCA loss function to find the top 100 forecasts from the same 10,000 iterations. We’ll set the minimum mean and minimum standard deviation of epsilon to 0, assuming a perfect match. We’ll colour by average cost and compare to the Soft L1 method.

Comparing Top 1% of Models Using Soft L1 and PDCA Loss

We can see that the PDCA Loss function has improved the distribution of fits above and below the actuals, a behaviour that we would like to carry forward because it helps to describe a high and low fit, as well as a bestfit. Lastly, let’s look at the parameter distributions for the top 1% of models to get an idea of the solution space using the PDCA Loss method.

Parameters of P1 model matches using PDCA Loss method

Summarize observations from these distributions.. They are mostly lognormal. Therefore, do we need to step into MCMC to get a proper PDCA distribution going forward? Keep these in mind and compare later on..

Summarize findings in these two charts and how we’re going to carry PDCA loss going forward in the rest of the evaluation / process as it provides a better probabilistic representation. Show distribution fits of the accepted parameters? Or do this at a later time? Or here just FYI??

Burn-In

The burn-in process is an initial convergence step that allows the MCMC to start from an equilibirum position, in this case the bestfit from the 10,000 Monte Carlo runs. While the burn-in process is often described as not needed or biased, in the case of PDCA we would like maintain its functionality to give us a baseline fit for evaluation of the PDCA loss function. In plain words, there is more uncertainty in production data that is hard to fit a curve to rather than data that is easy to fit. In this way, we can honour that higher level of uncertainty when needed.

Metropolis-Hastings Algorithm

With the bestfit defined as the starting point, we can start sampling through the solution space to carry out the MCMC process. The Metropolis-Hastings algorithm is a method for obtaining a sequence of random samples from multi-dimensional distributions like we have in PDCA. Simply defined, the Metropolis-Hastings method starts at an arbitrary multi-dimensional point (the start of the Markov Chain) and randomly chooses the next multi-dimensional point to evaluate given jumping distributions for each dimension. In our case the inputs parameters for the MH ARPs equation make up each dimension. For each iteration, the candidate sample is compared to the previous sample using an acceptance ratio. If this ratio is greater than 1, the jump is automatically accepted as being more probable than the previous sample. If the ratio is less than 1, the point is less probable than the previous sample, but can still be accepted. Using a random a random probability, sometimes the move will be accepted and sometimes it will be rejected, with acceptance ratios closer to 1 being given a higher probability of acceptance. Thus, we will tend to stay in high density regions and define the desired distribution after sufficient iterations.

The jumping distributions we will use are defined as follows:

Use some formulas to describe the N(i-1,SD) of each jumping distribution

Applying these jumping distributions for 10,000 iterations on our sample well, we can follow the chain of accepted or rejected iterations as the posterior distribution is defined.

Accepted input parameters using MCMC method

The resulting accepted distributions are dependent on the jumping distributions, which need to be tuned appropriately, but we can see the posterior distributions take shape. With more samples than the Monte Carlo illustration, each distribution is well defined, following normal or lognormal distributions.

Add note about staying in the proper solution space (prior distributions) through out entire MCMC process - this is crucial to have the equations make physical sense

Results

Let’s compare the top Monte Carlo matches to the MCMC matches to see how the accepted solution space matches the oil rate behaviour.

Comparing Top 1% Monte Carlo Matches to MCMC Matches

We can see that the matches have a tighter spread around the actuals, returning more forecasts within the acceptable solution space. The MCMC and Metropolis-Hastings algorithm method has achieved the desired result, returning thousands of acceptable models that define the uncertainty in forecasting the sample well. Using the same methods derived earlier in the markdown, let’s compile the accepted forecasts into P10, Mean, P50, and P90 durves to summarize the range of outcomes.

Finally, let’s summarize the EUR distribution and the EURs represented by the probabilistic forecasts.